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Abstract 

The presence of long-ranged correlations in a fluid undergoing uniform shear flow is investi- 
gated. An exact relation between the density autocorrelation function and the density-mometum 
correlation function implies that the former must decay more rapidly than 1/r, in contrast to pre- 
dictions of simple mode coupling theory. Analytic and numerical evaluation of a non-perturbative 
mode-coupling model confirms a crossover from 1/r behavior at "small" r to a stronger asymptotic 



power-law decay. The characteristic length scale is £ ~ ^J\$/a where Ao is the sound damping 
constant and a is the shear rate. 

PACS numbers: 05.20.Jj, 05.65. +b, 82.70.Dd 
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I. INTRODUCTION 



Long range spatial correlations (algebraic decay) in simple classical fluids at equilibrium 
occur only near a thermodynamic critical point, i.e. for finely tuned values of the thermody- 
namic parameters. On the other hand, such long range correlations appear generically for a 
wide class of nonequilibrium states |l[ . The predictions of this phenomenon have been made 
in a number of contexts, including self-organized criticality 0, linear response ||, kinetic 
theory and stochastic hydrodynamics ||. The simplest and most direct approach is that 
of linear response where the nonequilibrium correlation functions are expanded about the 
equilibrium state to first order in a nonequilibrium control parameter (typically a spatial 
gradient). The algebraic decay is then seen to result from spontaneous excitations of hydro- 
dynamic modes induced by the coupling of the control parameter to an associated flux. A 
more intuitive analysis of this effect follows from an extension of fluctuating hydrodynamics 
to nonequilibrium states. To linear order in the control parameter the same algebraic decay 
is found, as expected. Such theoretical studies for the density autocorrelation function in a 
fluid subject to a temperature gradient have received detailed experimental confirmation in 
recent years H] . The shortest length scale is set by the intermolecular force range, while the 
experimental verification is on macroscopic system size scales. Since the decay is algebraic 
there would appear to be no other length scale involved. However, we argue here that there 
is an additional macroscopic scale set by the parameters of the nonequilibrium state such 
that the true asymptotic decay is faster than that predicted by simple perturbative studies 
near equilibrium. The analysis here is limited to a single nonequilibrium state, that of uni- 
form shear flow, but the qualitative features are expected to extend to other nonequilibrium 
states as well. 

Uniform shear flow (USF) is characterized by a constant average density and temperature, 
and an average velocity flow field given by if ("r > ) = *o* • ~r* where the shear rate tensor 
is traceless ( *a* = ayx in a Cartesian frame of reference). This nonequilibrium state 
has a single scalar control parameter a and has been the subject of numerous theoretical 
investigations aimed at understanding transport and fluctuations in a model nonequilibrium 
state [0, |8|, |9|, [10], [TT| . All of these share in common the assumption that at large length and 
time scales, the dynamics of fluctuations in a simple fluid are dominated by the contribution 
of the hydrodynamic modes which decay much more slowly than do the neglected, kinetic 
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modes. The result is that the decay of thermal fluctuations in the hydrodynamic fields, i.e. 
the density, momentum and energy fields, is governed at large length and time scales by 
equations formally identical to the phenomenological Navier-Stokes equations. Here, "large 
length and time scales" means scales large compared to the mean free path and mean free 
time which is the usual domain of validity of hydrodynamics. Correlations between the 
values of thermally generated fluctuations in the fields at two different space-time points 
can be modeled by supplementing these equations with random forces which represent the 
interaction of this restricted set of variables with the neglected degrees of freedom to give 
a Langevin model. In equilibrium, the result of the Navier-Stokes-Langevin model is that 
the equal-time correlation functions for two hydrodynamic fields x(T*,t) and y(T*,t) are 
simply proportional to delta-functions in the spatial arguments 

C xy ( r>, !»') = (Sx ( T>, t) 6y ( T*', t)) -> A6 ( r> - r>') (1) 

where 5x = x — (x), and the amplitude A is a corresponding thermodynamic response 
function. This result simply confirms that fluctuations at different points in space are 
uncorrelated when speaking of hydrodynamic length scales (i.e., neglecting correlations on 
the scale of the force range). In contrast, these correlation functions for the nonequilibrium 
state have a new long range component, which to first order in the shear rate is of the form 

C xy (T>, r») - A6 (T* - !*') + aB \ (2) 

The amplitude B is again a thermodynamic response function. The physical difference 
between equilibrium and nonequilibrium is due to the way in which small hydrodynamic 
fluctuations decay. At equilibrium this occurs locally due to viscous and thermal damping, 
whereas in shear flow it is convected as well and spreads out over a length scale that varies 
as the speed of convection times the time scale for viscous and thermal damping M, . 
With the exception of reference , the form @ was generally obtained using a pertur- 



bative treatment in Fourier representation assuming that the shear rate a (with dimensions 
of frequency) is smaller than all other hydrodynamic frequencies, ck and Xok 2 , where c is a 
propagation velocity, Ao is a transport coefficient, and k is the wavevector. For fixed shear 
rate, this therefore sets an upper bound on the range of separations in real space for which 
the results apply. In this context, the true asymptotic behavior of the correlation functions 
remains unclear. The analysis of reference |T(J has the potential to resolve this question 



since it is non-perturbative and retains the two dominant effects of large shear rate: secular 
effects ~ at associated with convection, and shear rates comparable to the hydrodynamic 
damping a ~ X^k 2 . The resulting correlation functions are found to have the form 

C xy ( r», T*') -> AS ( T> - T") + aB 1 F ((T* - ~f*') /I) . (3) 

The function F(~r > /1) depends on a new nonequilibrium correlation length, I = y/\ Q /a . 
For r « I the result (^[) is recovered. However, the asymptotic form for r » I was not 



explored in any detail in reference [10] 



Here, we attempt to clarify the asymptotic behavior of the static correlation functions in 
a sheared fluid in two ways. First, it is shown that the continuity equation and stationarity 
place exact constraints on the decay of the density autocorrelation function such that it must 
be faster than r _1 for large r. This result is independent of any model for evaluating the 
correlation function. Next, we reconsider F(~f*/l) in (|3|) from the results of reference flO 



and show the actual asymptotic behavior is r~ n / 3 . The crossover between the r _1 behavior 
at short length-scales and the stronger r -11 / 3 decay at large separations is illustrated by 
numerical evaluation of the general result. 

II. EXACT BOUNDS ON THE RATE OF DECAY 

Consider N atoms with positions and momenta denoted ifj and ~p*i respectively and de- 
noted collectively as T = {ifi, ~p*i} i=1 . The atoms interact via a central two-body potential 
4> (r) and are confined to a volume V such that the average density is n. The potential is 
assumed to be repulsive at short distances and to diverge as r — > and to have a finite force 
range. Uniform shear flow results from the application of Lees-Edwards boundary conditions 



| 12|1 consisting of periodic boundaries in all directions except that of the gradient (here, the 
y-direction). If a particle exits the volume in the positive y-direction (say, at y — L/2), it is 
re-entered at the opposite side of the volume ( y = —L/2) with its velocity in the direction 
of flow shifted as v x — > v x — aL. Taken together, these constitute periodic boundaries in 
the local rest frame. The temperature will generally increase due to viscous heating but we 
will follow standard practice and assume the presence of a thermostat which counteracts the 
viscous heating so that the temperature is also constant in time leading to a steady-state 
13| . The dynamics thus described possess several important symmetries that will be used 
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below. First, they are invariant with respect to parity of the positions and momenta so that 
the equations of motion and boundary conditions are the same if written in terms of the 
variables = —ifi and ~p*i — — ~pV Second, they are invariant under mirror reflection 
about the z-axis (but not the sort/ axes since the boundary condition couples the x-velocity 
and the ^-coordinate). Third, since USF is a steady state, statistical properties are invari- 
ant with respect to a change in the origin from which time is measured (time translational 
invariance). Fourth, USF is translationally invariant in the local rest-frame as well as in a 
mixed frame consisting of the laboratory positions and the velocities measured relative to 



the flow [10 



The microscopic local density and momentum fields are defined respectively as 

N 



Vv^rtoro)) = ^(^-TiW) 

i=i 

N 

? p ("r*; T(t; r )) = ^ ^(t)6 (T* - ^(t)) (4) 



where ~p*i = — m i < li > ■ ifi is the momentum defined relative to the local flow field 
and r(t;To) is the point in phase space the system would reach after evolving from the 
initial point r for a time t. Note that these fields are related by the microscopic continuity 
equation 

d N d 

s *.(f ! r(t ! ro)) = 5:W)- 5f ^(f-W)) 

i=l 

= - V • (? P (t> ; r(* ; r )) + *o* • rVn (t*; r(t ; r ))) . (5) 

The correlations functions are defined as 

Cab r"; t, £') = / dT p {T ) 5^ a (T*; T(t; r )) (T* 7 ; r(f; r )) (6) 



where subscripts a and b label the specific field considered, p (r ) is the distribution of phase 
variables at initial time t = 0, and 5ip a = ip a — (if) a )- F° r a steady-state the time dependence 
occurs through t — if due to time translational invariance. Because of the modified spatial 
translational invariance, it is also possible to show that the correlation functions depend 
only on the relative separation. Consequently the correlation functions can be written 

C ab ( r> , r"; t, t') = C ab (ft - r*'; t - f) 
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Choosing T >/ = if and t' = t, Eq.(§|) gives the relationship 

^C nn ( T*) = = -V ■ (V pn (r>) - ~6 np (r>) + V • rC nn (r>)) (7) 

where for notational simplicity C n „ (T*) = C n „(T*;0). This result has been derived pre- 
viously in a different context [[Kj. It simplifies further using the translational, parity, and 
reflection invariance noted above to give the final form of interest here 

= -V ■ (2~tf pn (T*) + *a* ■ T>C nn ( r>)) . (8) 

This is an exact result that follows directly from stationarity and conservation of mass. 
Integrating (|j) over a spherical volume bounded by shells at r = + and r = R, and making 
use of Gauss's theorem gives the relationship of ( r ) to C nn ( t*') 

-2 y ^ p „ (i2r) • fdr = ^ f • • ~~r*C nn (Rr) dr (9) 

The notation cZr indicates a surface integral over the unit sphere and use has been made of 
the fact that the correlation functions evaluated at the origin vanish since the potential will 
not allow atoms to occupy the same spatial position (the singular contribution to C nn oc 5 (r) 
is excluded from the integration volume). 

Equation (|9]) is the main result of this section. To put this in context, consider an 
expansion of C nn (~r*) to first order in the shear rate with the form 

C nn (T>) - Cjg (r) + f ■ V • fCffi (r) + o(a 2 ). (10) 

Inserting this into eq.(H) gives 

-2 J ~6 pn (Rr) -rdr^ a^RC^ (R) + ( a 2 ) , (11) 

showing that if the density-momentum correlation function decays for large separations, 
as it must, then the first-order correction to the density- density correlation function must 
decay faster than 1/R . This result is independent of the dimensionality of the system. 
Furthermore, multiplying eq.(^) by R D_1 and integrating gives, in D-dimensions, 

J Jtt 

= a- J r D CW(r)dr + o(a"). (12) 



The quantity on the left is 




which vanishes in equilibrium but can be finite for USF due to velocity correlations. This 
in turn implies that the first order correction to the shear rate decays faster than l/r D+1 . 

These results can be given a somewhat more general interpretation. Any function of 
a vector in terms of the spherical harmonics in order to separate the dependence of the 
function on the direction and magnitude of its argument. The density-density correlation 
function becomes 

oo L 

Cnn ( T>) = E C nn (0 Y LM {?) . (14) 

L=0 M=—L 

with 

(r) = J d?Yt M (f) C nn ( r>) . (15) 

Because of the parity symmetry of USF, it is easy to show that only coefficients with even 
values of L are nonzero while the inversion symmetry about the z-axis implies that only 
even values of M contribute. Then, noting that 



Eq.(^|) becomes 




The conclusions drawn about the first-order correction C n n (r) are seen to be exact state- 
ments also about lmC^ n (r) , valid to all orders in the shear rate: namely, that lmC^ n (r) 
must decay faster than 1/r in any number of dimensions and faster than \/r D+1 in D - 
dimensions if the spatial integral of the radial part of the density-momentum correlation 
function is finite. This corresponds precisely to the quantity for which the first-order results 
(0) predicted a 1/r decay. Consequently, that result cannot be correct for sufficiently large 
r . 
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III. APPROXIMATE EVALUATION OF C nn (T) 



The algebraic decays of both static and dynamic correlation functions observed in simple 
fluids can be derived by means of fluctuating hydrodynamics. In this model, the exact 
conservation laws for the local density, momentum and energy density fields, p, J? and u 
respectively, 

|p> + S?.<P> = 
JU+V'-Tf = (18) 

are approximated by taking the pressure tensor, P , and heat flux vector, if, to be a sum of 
two terms: the usual Navier-Stokes functionals of the local fields, and a random component 
which is delta-function correlated in space and time. The amplitudes of these correlations 
are related to the forms of the deterministic parts of the fluxes |14|]- Rewriting these in terms 
of deviations from the macroscopic state then gives a description of fluctuations about this 
state. The details for USF have been discussed in detail elsewhere and only the results 



are quoted for the purposes here. Defining the Fourier transform of the density-density 
correlation function as 

C nn (is^j = J d~r* exp (f£ ■ T^j C nn (T*) (19) 
the result obtained for it is 

Cnn (y-, a) = k B T plxT (l + 7" 1 A„„ ) (20) 



A„„ (¥j = J°° ds kk ^} s) S) exp (- J ds' k 2 (-.')] ■ (21) 

Here xt is the isothermal compressibility, cq is the equilibrium speed of sound, Ao is the equi- 
librium sound-damping constant, 7 = c p /cy is the ratio of specific heats at constant pressure 



and volume, and k (t) = (k x ,k y — tk x ,k z ). The characteristic length scale is / = a/Ao/o 
. The first term in eq. (|2~0D is the equilibrium contribution which, in the small wavevector 
approximation used here, is a constant. The second term represents the nonequilibrium cor- 
rection which is derived assuming that the shear rate and rate of dissipation are significantly 
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less than the sound frequency, a, \ k 2 <C c$k where Co is the speed of sound. However, no 
restriction on the value of a/Xok 2 is imposed. The inverse transform of eq. (^) is the real 
space result quoted in ([3]) above. 

In order to evaluate the behavior in real space, it is useful to introduce the expansion of 
C nn [ k ) in spherical harmonics 



oo L 



C nn (*")=££ C<nf (*) Y LU (k) 
L=0 M=-L 

and to recall the relation between the coefficients ([L5|) to those in ( 22 ) 



(r) 



2txH 



2aL 



k 2 dk 3L (kr)c::(k). 



(22) 



(23) 



This gives directly 

C (r) = (fc fl T pgx T ) v^r^ ( r/0 SloSmo + (Mp 2 ^" 1 ) ^_r s A^ (r/Z) (24) 
with 



Kn (p = r/l) 



dk' j L (kp)Yl M (k) / ds kk ^ s) 3) 



cxp 



ds' fc 2 f-s r 



x 



dfc jl (fc) y^ M ( k 



exp (~ f Q ds ' fc 2 (-sV)) . 



(25) 



/o k 3 (-sp 2 ) 

The small shear rate limit (|2|) is obtained by noting that I 2 is inversely proportional to the 
shear rate and expanding to leading order in a 



r 3 Ai M (r/l) - I 



- I dk' j L (k)Y* M [k)^ 



a 

\ r 



i$L2 (Sm2 - S M -2) 



7r 2n 




(26) 



4 V 15 

The context here shows that this result applys only in the limit r/l — > which is to say 
that, for fixed shear rate, this is a small r result and does not represent the true long-range 
behavior of the correlation function. 

The actual asymptotic behavior for large p = r/l is obtained in the Appendix where it is 
found that 



AS (P) 



2/3 / 

()j svrr Qj ) p-"/3 + o (p-1 3 / 3 ) 

6.5p- n / 3 . 



(27) 





FIG. 1: In (A^) as a function of ln(r//) as determined by the numerical calculation (circles) and 
the asymptotic result given in eq,([27|) (line). 



All other components decay even more rapidly. In particular the true asymptotic behavior 
of ImA??, (x) is found to be 



ImA- (p) - 9 ^Q 1/6 V5rrT + O (p 



-19/3 



) 



13.2 p 



-17/3 



(28) 



This is consistent with the exact result (\L7\) that this component must decay more slowly 
than r -4 in three dimensions. 

In order to probe the asymptotic behavior in more detail, we have performed a numeri- 
cal evaluation of eq . (p5|) by means of multi-dimensional Monte Carlo integration using the 



VEGAS algorithm [16], [17], [18). Rather than directly evaluating eq. (|25|) , it was found to be 
more efficient to separate out the short-ranged 1/r behavior by rewriting this as 

'27T 7T 



. LM 



(P) 



i$2l ($m2 - 5 m _ 2 ) 



+ Jdt jl (kp) Yl M (k 



x / ds 
'o 



k ( s) k x ky 



cxp 



ds' k (-0 



(29) 



k 3 (-s) k A j x ( , 

The number of samples used in the performing the integrals was adjusted so that the internal 
estimate of the error in the evaluations was always less than 5% of the calculated values. 
For small p, the errors were substantially less while the limit was occasionally reached as 
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FIG. 2: In (JmA^j as a function of ln(r//) as determined by the numerical calculation (circles) 
and the asymptotic results given in eq.(p8|) (full line) and eq,(p6[) (dotted line). 



p was increased. Figure 1 shows the spherically averaged value, A°° (p), as a function of 
p together with the asymptotic power law of —11/3 and the two are seen to be consistent. 
Figure 2 shows the numerical calculation A^ (p) in comparison with the small- and large-r 



limits, eqs. (|26|) and ( p8|) respectively. Again, crossover between the limiting forms is clearly 
identified. 



IV. DISCUSSION 

The prediction of 1/r decays in the density autocorrelation function in shear flow is in 
violation of an exact bound coming from elementary considerations of statistical mechanics 
and the properties of the USF steady state. Detailed analysis of a less restrictive solution of 
the Navier-Stokes-Langevin model confirms that the 1/r behavior is actually valid only for 
r < a/Ao/o and that this crosses over to a stronger power law decay at large distances. It 
turns out for this model that the density-energy correlation function as well as energy and 
longitudinal-velocity autocorrelation functions share the same spatial dependence given by 
(p6|). Therefore the calculations given here apply to them as well. Interestingly, the most 
long-ranged correlation function, based on calculations analogous to that illustrated in the 
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appendix, is for one of the transverse- velocity autocorrelations which, in the notation of ref 
PTCj] , has the form 

C 44 (T>) = k B T [1 + A U (T>/1')] (30) 
with V = \j2v$/a where vq is the shear viscosity and 

A™( P )^lQV\r(^p-^ + 0(p-^). (31) 

The analysis here has implicitly assumed stability of the USF state. In fact USF is 
unstable to sufficiently long wavelength perturbations [11, [15|]. The critical wavelength for 
stability scales approximately as Vt/ol for small a where Vt is the thermal velocity. Therefore 
in order to see the crossover phenomenon discussed here there must be conditions such that 
/ << vt/cl. This requires aXo/v^ << 1 which can be accomplished by small shear rates 
and high temperatures. The predicted behavior should be accessible via molecular dynamics 
simulation for a sufficiently large system. 

The qualitative feature of a nonequilibrium length scale should be more general than the 
special case of USF, and applicable to other nonequilibrium states as well. For example, 
a steady state with uniform temperature gradient has a characteristic frequency u^V In T. 
Setting this equal to a hydrodynamic damping gives the length scale I = y/Xo/ (^tV InT). It 
is expected that the asymptotic decay for r » I will be different from that of perturbative 
mode coupling theory currently in the literature for reasons similar to those given here for 
USF. 
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APPENDIX A: ASYMPTOTIC BEHAVIOUR OF THE CORRELATION FUNC- 
TIONS 

We begin by writing 

~ f o f°° k k -+- tk 2 

A nn (t)= dt- + %K exp {-k 2 p(t)) (Al) 

" ' Jo (l + 2k x k y t + k 2 x t 2 j 
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where (3(t) = t + k x k y t 2 + \k 2 t z is independent of the magnitude of the wavevector. The 
coefficients of the expansion in spherical harmonics in real space are then 



LM 



(p) 



dkYZ M (k)J dkk 2 3l (kp) J dt 



k x ky -\- tk x 



(l + 2k X kyt + kit 2 ) 



3/2 



exp (-k 2 (3(t)) 



[ 2 ] -p- 3 



2 L + 2 V (2£±3) 



X 



k x ky -\- tk x 



(i+2k x k y t+k 2 x t 2 y /2 



L+3 



M 



L + 3 2L + 3 



2 ' 2 ' Af3(t) 



A2) 



where M(a, b, z) is the confluent hypergeometric function. Now, changing variables to t 



2/3 



(3) 



-1/3 



y gives 



LM 



(p) 



X 



2 L+2 r (2L^3) 

/ 

./0 



r-^ 3 I dkYl M (k){kl) 

-2/3 



-1/6 



1 + k x ky ( /c^ 



2/P 



-2/3 



l + 2k x k y (kiy 2/3 yp-^+(kl ) ) 



-1/3 \ 3 / 2 

~ 5 - 4 /3 



1/ 2 P" 



X 



Vt(s/) 



3 \ 2 



M 



L + 3 2L + 3 



with 



7(y) = 3 + kxky ( k 



-2/3 



2/P 



- 2 / 3 + u 



4 7 (y) 

1/3 



1/ 2 P" 4/3 . 



(A3) 



(A4) 



For the 00 component, we use M (a, a, z) = exp(z) to get 



23' 



- 11 / 3 dk k 2 



-1/6 



cfy (3y 3 ) 2 exp 



3y J 



+ O (p~ 13 / 3 ) (A5) 



= 2§3-5 5vrr ^ p- 11 /' 3 + O (p- 13/3 ) . 



For the 22 component, it is more convienient to go back to the beginning and to integrate 
by parts 

A£f (p) = f dk Y* LM (k) J™ k 2 dkj L (kp) 



x 



= P 



J dt -j t (^ + 2 k x k y t + k 2 x t 2 ^j 
dkYl M (tjJ k 2 dkj L (k)-V^ 



exp (-k 2 p(t)) 



-p 



Jo 



dt (l + 2k x k y t + k 2 x t 2S ) 



1/2 / p 



Pit) 



2 L + 2 Y (L + |) 

L + 5 2L + 3 



dk Y* LM (k) 



M 



2 ' A(3(t) 



A6) 
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The first term vanishes for L > 1 so 



x * (l + 2k T k,t + J&) V2 2 exp (-J^) (AT) 



and making the same change of variables as above gives 



ImA^(p) = -^P~ U/3 J dkl™Y* 2 (t) (%l 

x J™ y-*dy (l + 2^p" 2 / 3 (hi) " 2/3 y + p~ 4 / 3 (kl) y 2 ) 



1/2 



x (i) %xp (-4) (A8) 

Now, only the odd terms in k x k y give nonvanishing contributions so we expand as 

(l + 2^p-V3 ( S ) ~ 2/3 y + p -V3 ( g ) ^ * ^ * ^ 

= ^V3fcA (S) _2/3 (% 3 - 38) exp (~^ 3 ) p^ 3 + ^Zyi% x k y (fc^exp (-jjy 

x ^8 (£ 2 ) (-918y 3 + 81y 6 + 2008) + (k x k y ^ 2 (28 04% 3 - 5022?/ 6 + 243y 9 - 40 088) 
+0(r~ 8/3 ) + even 

The ^/-integral of the first term vanishes and the next term gives 



x / ( //,- I ui^ ( k ) ( k 2 x ) " 1/6 k x k y (kA " 2 ( 3 ( fc 2 ) - 2 feC ' 



^3!V302f «T g) p- 17 / 3 (A10) 
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